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ABSTRACT 

We argue for implementing star formation on a viscous timescale in hydrodynam- 
ical simulations of disk galaxy formation and evolution. Modelling two-dimensional 
isolated disk galaxies with the Bhatnagar-Gross-Krook (BGK) hydrocode, we verify 
the analytic claim of various authors that if the characteristic timescale for star forma- 
tion is equal to the viscous timescale in disks, the resulting stellar profile is exponential 
on several scale lengths whatever the initial gas and dark matter profile. This casts 
new light on both numerical and semi-analytical disk formation simulations which 
either (a) commence star formation in an already exponential gaseous disk, (b) begin 
a disk simulation with conditions known to lead to an exponential, i.e. the collapse 
of a spherically symmetric nearly uniform sphere of gas in solid body rotation under 
the assumption of specific angular momentum conservation, or (c) in simulations per- 
formed in a hierarchical context, tune their feedback processes to delay disk formation 
until the dark matter halos are slowly evolving and without much substructure so 
that the gas has the chance to collapse under conditions known to give exponentials. 
In such models, star formation follows a Schmidt-like law, which for lack of a suit- 
able timescale, resorts to an efficiency parameter. With star formation prescribed on 
a viscous timescale however, we find gas and star fractions after 12 Gyr that are 
consistent with observations without having to invoke any "fudge factor" for star for- 
mation. Our results strongly suggest that despite our gap in understanding the exact 
link between star formation and viscosity, the viscous timescale is indeed the natural 
timescale for star formation. 
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1 INTRODUCTION 

The origin of the exponential radial light profiles in galaxy 
disks is controversial. The in situ scenario argues for stars 
forming in an already exponential gaseous disk (Gunn 1982; 
van der Kruit 1987; Dalcanton, Spergel & Summers 1997; 
Mo, Mao & White 1998) . The viscous scenario, on the other 
hand, contends that regardless of the initial gas density pro- 
file in the disk prior to star formation, the stars will acquire 
an exponential profile, provided star formation proceeds on a 
viscous timescale (Silk & Norman 1981; Lin & Pringle 1987; 
Clarke 1989; Yoshii & Sommer-Larsen 1989; Saio & Yoshii 
1990; Olivier, Blumenthal & Primack 1991; Silk 2001; Fer- 
guson & Clarke 2001). 

The in situ scenario defers the problem of the origin 
of exponential stellar light profiles to the origin of the ex- 
ponential gas distribution. The explanation for such a gas 
distribution stems from a remark made by Mestel (1963) 



on the similarity of the specific angular momentum distri- 
bution observed in the disks of spirals to that of a uniform 
sphere in solid body rotation. If specific angular momentum 
is conserved during the collapse of such a gaseous sphere, 
the resulting disk naturally retains the sphere's original an- 
gular momentum distribution. Considering the case where 
this uniform sphere of gas collapses within an isothermal 
dark halo, Gunn (1982) found that the gas surface density 
distribution in the final disk would be approximately expo- 
nential over three scale lengths. 

This simple model is the basis of recent explanations 
of a host of disk properties. By extending Gunn's scenario 
to allow the dark halo to adiabatically respond to the ac- 
cumulation of baryonic mass in its center, Dalcanton et 
al, (1997) explained not only the structural properties of 
present-day disks {e.g. their exponential stellar profiles) but 
also their dynamical properties {e.g. their rotation curves, 
the TuUy- Fisher relation). Using the Press-Schechter pre- 
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scription to compute the abundance of virialized halos (Press 
& Schechter 1974), combined with the prediction from nu- 
merical simulations of both the halo dark matter density 
profiles (Hernquist, 1990; Navarro, Frenk & White 1996) 
and spin parameters (Barnes & Efstathiou 1987; Warren et 
al. 1992), these authors also predicted the surface bright- 
ness and scale length distributions of present-day disks. Mo, 
Mao & White (1998), using a similar modelling, extended 
this study to predict the evolution of disk properties within 
various cosmological models. 

In addition to the simple initial conditions set out by 
Mestel and Gunn, models built upon the m situ star for- 
mation scenario require a recipe for converting gas to stars. 
In the case of Dalcanton et al. (1997), Mo et al. (1998), and 
more generally in all semi-analytic models as well as numeri- 
cal simulations, a Schmidt law, with an efficiency parameter 
is adopted (see e.g. Kauffmann, White & Guiderdoni 1993; 
Cole et al. 1994; Somerville & Primack 1999; Katz 1992; Mi- 
hos & Hernquist 1994; Steinmetz & Miiller 1995; Gerritsen 
& Icke 1997; Springel 1999). Of course, to imprint the expo- 
nential profile on the stars, such a prescription must ignore 
any process that could significantly redistribute the gas mass 
and angular momentum in the disk, such as gaseous outflows 
or viscous flows induced by a turbulent ISM or spiral arms. 

Numerous hydrodynamical disk formation simulations 
have shed light on the in situ theory for the origin of expo- 
nential disk profiles. The earliest were N-body plus SPH sim- 
ulations, first without star formation (Katz & Gunn 1991) 
and then with stars (Katz 1992). These simulations followed 
the collapse of an isolated, uniform density perturbation 
in solid body rotation embedded in an expanding universe. 
External tidal fields were neglected. Since some small-scale 
power was added, dark matter and gas clumped prior to the 
gas dissipating its energy and falling to form a disk. Obvi- 
ously, the merging of these subclumps involved angular mo- 
mentum transfer, but apparently the deviation from Mestel 
and Gunn's idealized initial conditions was minor and the fi- 
nal gaseous and stellar profiles were well-fit by exponentials. 
All other simulations starting from similar initial conditions 
have led to the same conclusion (Vedel, Hellsten & Sommer- 
Larsen 1994; Steinmetz & Miiller 1995; Contardo, Steinmetz 
& Fritz- von Alvensleben 1998). 

In the current paradigm of hierarchical galaxy forma- 
tion however, one needs to follow the formation of disk galax- 
ies in the context of a more realistic tidal field, for it seems 
implausible in such a paradigm that all gas destined to set- 
tle into a disk originated in a rigidly rotating spherical pro- 
tocloud which conserved its specific angular momentum not 
only during the process of protogalactic collapse and settling 
into centrifugal equilibrium but also throughout the star for- 
mation process. Unsurprisingly, hydrodynamical simulations 
of hierarchical galaxy formation illustrate how nonspherical 
subclumps of gas and dark matter torque eachother while 
merging. In fact, in such simulations of only gas and dark 
matter, there is so much transfer of angular momentum dur- 
ing merging that the disks which form have angular mo- 
menta two orders of magnitude lower than what is observed 
(Navarro & Steinmetz 1997). 

To remedy this, attempts have been made to delay the 
cooling and subsequent collapse of gas in dark halos, so that 
disks could form late in cosmological simulations when dark 
matter halos are evolving more slowly and have much less 



substructure (Weil, Eke & Efstathiou 1998). In these situ- 
ations, the idealized initial conditions envisioned by Gunn 
(1982) may be an accurate starting point for disk forma- 
tion. First attempts at physically motivating such a delay 
by including feedback mechanisms instead of just switch- 
ing on cooling in a simulation at a late epoch have been 
undertaken by Sommer-Larsen, Gelato & Vedel (1999). In 
their most successful cases, the gas surface densities of their 
disks are poorly represented by exponentials (their Figure 
4). However, resolution in their simulation (between 3303 
and 5017 gas particles in a disk) is such that determining 
a disk scale length is already a difficult task. Hence it may 
be resolution that prevents these profiles from looking more 
exponential. At any rate, on the evidence of current numer- 
ical simulations, we find it hard to be convinced that all or 
even most gaseous disks forming in a hierarchical context 
will acquire an exponential profile. 

Rather than trying to fine tune the simulations to ob- 
tain idealized initial conditions, it seems more appealing to 
look for a mechanism that could drive the stellar profile to- 
wards an exponential, regardless of the gas profile it starts 
from. Bearing in mind that star formation and the accompa- 
nying complexity of stellar feedback are poorly understood 
and would be hard to simulate (for obvious resolution is- 
sues) even if they were understood, the best that one can 
do is to model processes which govern star formation on a 
large scale. With this intention, as an alternative to Katz's 
(1992) star formation recipe based on the physics of Jeans 
instability, we propose to capture the physics that drives star 
formation by modelling the interstellar medium as a turbu- 
lent fluid of colliding, fragmenting and star-forming clouds 
with energy to drive the turbulence supplied from super- 
novae, stellar winds, and ionizing radiation from newly born 
stars. Such turbulence would drive viscous flows. 

Therefore, as a start, in this paper, we use numerical hy- 
drodynamical simulations to revisit the hypothesis that star 
formation occuring on a viscous timescale in galaxy disks 
yields exponential stellar profiles. Lin & Pringle (1987) ex- 
plored this question by numerically integrating ID differ- 
ential equations, and Yoshii & Sommer-Larsen (1989) con- 
firmed their result by deriving analytical solutions for the 
star forming viscous disk. Saio & Yoshii (1990) extended 
these works by including self-gravity. The numerical hydro- 
dynamical approach finds its justification in the fact that, 
in the future, we want to extend these investigations beyond 
ID calculations, and with physics that is difficult to include 
in an analytic approach {e.g. stellar mass loss, supernovae 
feedback) . 

The outline of the paper is as follows; after a discussion 
of turbulence as a plausible source of viscosity in galaxy disks 
(section ^) and a description of our numerical approach (sec- 
tion ^), we commence with benchmarking our code to the 
viscous evolution of a pure gaseous disk without star forma- 
tion (section ^). Then we proceed to a detailed examination 
of experiments of star formation in viscous galaxy disks with 
different initial gas profiles and dark halo rotation curves 
(section P). Finally, we discuss our results (section 
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2 PHYSICAL MOTIVATION FOR VISCOSITY 

A general consensus among the astrophysical community is 
that the true kinematic viscosity of the interstellar medium, 
v, is much too low to alter the galactic disk structure on 
a timescale comparable to the age of the universe. It has 
been argued by various authors (Lynden-Bell and Pringle 
1974; Silk and Norman 1981; Lin and Pringle 1987; Yoshii 
and Sommer-Larsen 1989; Clarke 1989; Olivier et al. 1991; 
Duschl et al. 2000; Silk 2001) however, that viscous tur- 
bulence could drive a non-negligible radial flow of gas in 
galactic disks. 

Turbulence is a consequence of fluid motion, therefore 
to be maintained it requires an external energy source E^xt- 
This energy which arises on large scales, with characteris- 
tic velocity U and length L, is progressively transferred to 
smaller and smaller scales by non-linear interactions, un- 
til the dissipation scale Lms is reached where molecular 
processes convert it into heat. It follows from dimensional 
arguments that Eext ~ U^/L. Similarly, one can obtain 
Ld^s ~ (v^E-J^ty^*- As a result, 

-^ = Re^/\Re=— (1) 

where Re is the Reynolds number. 

Although star formation and/or magnetic fields cer- 
tainly participate in feeding turbulence, we may suppose 
that turbulence is driven by differential rotation, (as argued 
by e.g. Richard and Zahn (1999) for accretion disks), for 
the following qualitative argument does not change. Typi- 
cal values for a galactic disk are L ~ lO^^m (radius of an 
orbit), U ~ lO^m s~^ (corresponding azimuthal velocity) 
and V ~ 10~^^m^ s~^ (we have assumed an average den- 
sity of 10^ at m ^, a temperature of lO^K and described the 
gas as elastic spheres of diameter 10~^°m). Combined to- 
gether, these numbers lead to Re ~ 10^^. This extremely 
high Reynolds number means that in order to fully resolve 
our problem numerically in 3D, we should have about lO'^^ 
grid points (see eq. |l|), whereas state-of-the-art super com- 
puters can handle 10^" at most. Therefore, the only way 
around this problem is to somehow model what happens on 
subgrid scales and numerically resolve the properties of the 
flow on large scales. 

As first noted by Boussinesq (1870), it seems that to 
first order and on macroscopic scales the net effect of tur- 
bulence is simply to turn the real fluid into a more viscous 
one. This can be schematically understood as follows: the 
Reynolds number of a turbulent fluid is so high that any 
perturbation in the smooth flow produces eddies which dissi- 
pate energy until the flow reaches a critical, lower Reynolds 
number Rcc, corresponding to the onset of turbulence. If 
more eddies are produced, the effective viscosity becomes 
too high, cutting off the production of eddies. On the other 
hand, if not enough eddies are produced, then the flow still 
exhibits a tendency to produce turbulence. Laboratory flow 
values for critical Reynolds numbers are on the order of 
Rcc ~ 10— 10'^ (see for example Richard and Zahn 1999), so 
that one can solve the Navier-Stokes problem with this crit- 
ical Reynolds number and hope that it is enough to obtain a 
good description of the properties of the flow on macroscopic 
scales. 

Our goal here is not to derive a mathematically rig- 
orous formulation of hydrodynamic turbulence in galactic 



disks but to suggest that there exist strong physical argu- 
ments in favor of describing the behaviour of the gas flow on 
large scales by simply solving a Navier-Stokes problem with 
a smaller Reynolds number, namely the critical Reynold's 
number, Rec ~ 10 - 10^. 



3 NUMERICAL APPROACH 

We perform numerical hydrodynamical simulations in 2D 
with the BGK (Bhatnagar-Gross-Krook) code (Prender- 
gast and Xu 1993, Slyz and Prendergast 1999). This is a 
scheme based on gas-kinetic considerations, meaning the 
code solves for the time evolution of a gas distribution func- 
tion, f{x,u,t), in phase space. Fluxes are computed from 
velocity moments of /. / evolves in time according to an ap- 
proximation to the coUisional Boltzmann equation, namely 
the BGK equation: = where is a time rate of 
change along the trajectory of a single particle moving freely 
in phase space under the action of a smoothly varying gravi- 
tational field and where models the effects of collisions. 
Here <; is a local equilibrium distribution, i.e. the Maxwell- 
Boltzmann distribution function, having the same mass, mo- 
mentum and energy densities as /, and r is a relaxation time, 
which can be as small as the mean time between collisions of 
a particle in the gas (Bhatnagar, Gross & Krook 1954). The 
underlying hypothesis of the BGK model is that there is a 
relaxation process in which collisions tend to reshuffle the 
particles defining the true distribution function, f{x, u, t) to 
g{x,u,t), while conserving the total mass, momentum and 
energy. This conservation requires that the moments of the 
BGK collision term must vanish, f's relaxation to the equi- 
librium state g is dissipative and hence is accompanied by 
an increase in entropy. It is nontrivial that the BGK scheme 
captures this physics. Many schemes have to invoke an "en- 
tropy fix" usually in the form of artificial viscosity to mimic 
an entropy increase and eliminate unphysical numerical phe- 
nomena such as rarefaction shocks. 

Together with the criteria of vanishing moments of the 
collision terrn, a numerical solution of the formal integral 
solution (eq. H below) of the BGK equation comprises the 
BGK scheme. Given initial conditions f = fo a.t t — to = 0, 
the solution to the BGK equation contains two terms: an 
integral over the past history of g and a term representing 
relaxation from an initial state, fo on a timescale t. 

f{xf,u,t)^- I 3(f ,u ,t )e~'*~* '/^dt -I- 
Jo 

e~''^^ fo{xf ~ ut,u,to) (2) 

Here x and u in the arguments of g are points along a 
gas particle's trajectory in phase space with the final con- 
ditions X = Xf and u = u at t = t. When fluxes are 
eventually computed at Xf, f{xf,u,t) is integrated over the 
full and continuous range of velocities, u, from — oo to +oo. 
Therefore fluxes computed from the true distribution func- 
tion f{xf,u,t) arise by a weighted integration of the equilib- 
rium distribution function g{x ,u ,t ) not just over one tra- 
jectory but over all possible trajectories (in 6-dimensional 
phase space) which arrive at (xj,t). A detailed description 
of the numerical procedure in the BGK scheme for the 2D 
Cartesian case is given in Xu (1998). 
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By evolving / through an equation which accounts for 
particle collisions, the fundamental mechanism for generat- 
ing dissipation in gas flow, the BGK scheme gives fluxes 
which carry both advective and dissipative terms. This 
means that for viscous problems no modification to the BGK 
flux expressions from their form for non-viscous problems is 
necessary. Furthermore there is no need for viscous source 
terms. More than this, unlike the earlier gas-kinetic scheme 
used in astrophysics, namely the beam scheme (Sanders and 
Prendergast 1974), which also carried dissipation in its flux 
expressions, the BGK scheme is dissipative in a realistic way. 
Where the beam scheme arbitrarily chooses the form of / 
at the beginning of each updating time step and evolves 
it through the coUisionless Boltzmann equation, the BGK 
scheme solves for the time evolution of / throughout an up- 
dating time step using a model of the coUisional Boltzmann 
equation. If the gas happened to be in disequilibrium at 
the beginning of the updating time step, the beam scheme 
cannot readjust the gas properly to equilibrium because 
it evolves / through the coUisionless Boltzmann equation 
which is simply missing the adjusting mechanism provided 
by collisions. 

Also by assuming instantaneous relaxation back to the 
chosen form for / at the beginning of an updating time 
step, the beam scheme endows the gas with a mean colli- 
sion time equivalent to the updating time step. Collisions in 
the BGK scheme on the other hand are active throughout 
the updating time step and for hydrodynamical applications 
the BGK scheme demands that the collision time is much 
smaller than the updating time step. Since dissipation pa- 
rameters are proportional to the collision time, r, e.g. the dy- 
namical viscosity r] — rp, and the heat conduction coefficient 
K = T ^2^ ^ (where p, k, m are the gas pressure, Boltzmann 
constant, and particle mass respectively and K = 
where 7 is the ratio of specific heats, — ), we easily see that 
an overestimation of the collision time will lead to a very 
diffusive scheme. 

The BGK code has been extensively tested on stan- 
dard ID and 2D test cases of discontinuous nonequilibrium 
flow with features ranging from cavitation (Roe test) to the 
collision of strong shocks (Woodward- CoUela test) (see Xu 
1998 for a review). At shocks and contact discontinuities 
the code behaves as well as the best high resolution codes 
which do not employ regridding in the neighborhood of a 
shock front and it gives better results at rarefaction waves. 
Because the BGK scheme naturally satisfies the entropy con- 
dition, it never requires an entropy fix when it encounters 
strong rarefaction waves (cf. the Sjogreen test). Tests of the 
BGK schemes ability to solve the Navier-Stokes equations 
in smooth fiow regions include the Kolmogorov and laminar 
boundary layer problem (Xu and Prendergast 1994). Tests of 
the BGK scheme in an external gravitational potential show 
the long-term stability of the scheme and its convergence to 
the equilibrium solution (Slyz & Prendergast 1999). 



3.1 Quantifying Viscosity in BGK 

As already mentioned in the previous section, on a micro- 
scopic level dissipative parameters are directly related to the 
collision time of the gas, r. Like other finite volume schemes, 
the BGK scheme follows the mass, momentum and energy 



densities within cells on an Eulerian grid. For the flux com- 
putation, the BGK scheme solves for / in the neighborhood 
of the boundary of each computational cell and for a short 
time (given by the usual CFL (Courant-Friedrichs-Lewy) 
condition). Hence at each time step and at each cell bound- 
ary, it is possible to solve the BGK equation with the colli- 
sion time T which is appropriate to the local fluid properties. 
For example, since gas kinetic theory gives a dependence of 
collision time on gas density p and temperature T, an expres- 
sion for r may be written as r = where Ci is a propor- 
tionality constant which is chosen according to the desired 
Reynolds number of the problem. (Ci — {U L / Re)[p'^ / y/\p) 
where U and L are characteristic velocity and length scales. 



A 



where T is the gas temperature, p is the gas pres- 



sure, and Re is the Reynolds number.) If we wish to do a 
problem with constant kinematic viscosity for the case where 
the Reynolds number is equal to the critical Reynolds num- 
ber, Rec, and an isothermal gas, then we set r to have a 
constant value, r = ■ Once r is known, we can write 

an expression for the real physical kinematic viscosity coef- 
ficient V = T-. 

p 

For problems in which the gas flow has discontinuities, 
the collision time r serves an additional purpose. If the grid 
is not fine enough to resolve a discontinuity, then artificial 
dissipation must be added to broaden the discontinuity so 
that it is at least one grid cell thick. Because viscosity and 
heat conductivity are proportional to r, the BGK scheme 
broadens shocks by enlarging r at the location of the discon- 
tinuities. The expression for the collision time in the BGK 
scheme therefore contains another term in addition to the 
real physical collision time. The second term is chosen in 
such a way that shocks in the fiow span at least one grid 
cell. In one-dimensional Cartesian coordinates this second 
term is: 



\Pl 



-Pt\ 



+ Pr) ' 



(3) 



The subscript I (r) denotes quantities interpolated from the 
left (right) side of the cell interface at which we compute the 
collision time r, and p is the gas pressure. 

This second term tunes the amount of artificial dissi- 
pation in the scheme. The notable difference between how 
the BGK scheme inputs artificial dissipation and how other 
schemes input artificial dissipation is that the BGK scheme 
puts it in exactly as if it were real dissipation corresponding 
to the numerically necessary value for r. We again emphasize 
that neither the artificial dissipation nor the physical dissi- 
pation enter the scheme as a source term. When the second 
term (eq. |^ in the collision time dominates, the order of the 
scheme is reduced and the true distribution function / is 
determined more by the initial distribution function /o than 
by the integral over g (eq. |^. For smooth fiow however, such 
as is encountered in the viscous evolution of a galactic disk 
in an axisymmetric gravitational potential, C2 may be set 
to zero and the BGK scheme evolves the gas with its true 
physical dissipation parameters, as long as the grid resolu- 
tion corresponding to the desired Reynolds number of the 
problem is achieved. 

Even when C2 is set to zero, in addition to the (ex- 
plicit) real physical viscosity in the scheme (controlled by 
the collision time) there is an inevitable numerical viscos- 
ity, dependent on the grid spacing: the coarser the grid, the 
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Figure 1. Early viscous evolution (first 2.28 Gyr) of a disk of 
pure gas initially in inviscid centrifugal equilibrium {vji = 0) in 
a fixed gravitational potential. Coding for the different curves is 
given in the bottom panel of the figure. 
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Figure 2. Same run as in fig. |^ but with the results for the evolu- 
tion shown up to 11.84 Gyr. The steady-state analytic predictions 
(asterices) are compared to the BGK numerical results (thin solid 
curves and empty circles). 



larger this numerical viscosity. The origin of this viscosity 
is that the scheme forgets everything except cell averages 
at the end of each time step. A consequence of this is that 
when one does a resolution study of the viscous evolution 
of a gaseous disk, the flow properties on a 207 X 207 grid 
quantitatively match those on the 85 X 85 grid only when 
one runs the problem on the finer grid with larger r. 



4 VISCOUS EVOLUTION OF A DISK OF 
PURE GAS 

Before testing the hypothesis that ~ t^isc gives exponen- 
tial stellar profiles, we check how closely the BGK hydrocode 
reproduces our analytic predictions for the steady state vis- 
cous evolution of a gas without star formation. The interest 
of performing such a computation is obvious: it allows us 
to measure the true viscosity of the code as well as its abil- 
ity to capture the Navier-Stokes solution for a compressible 
viscous fluid evolving in a fixed gravitational field for the 
case where the shear fiow is highly supersonic (Mach ~ 20) 
throughout most of the disk. Earlier tests of BGK's ability 
to give Navier-Stokes solutions were only for subsonic shear 
flows (Mach ~ .5 for the laminar boundary layer problem; 
Mach « .8 for the decaying sinuisoidal wave problem) (Xu 
and Prendergast 1994). 

We initialize the gas in a fixed dark halo with gravita- 



tional potential: 

Hr) = Jvo'ln(ro' + r') 



(4) 



where uo = 220 km/s and ro — 2.56 kpc. We consider two 
cases: the case where the gas is initially in inviscid centrifu- 
gal equilibrium (vh = 0), and the case where the gas is in 
centrifugal equilibrium but the initial radial velocity profile 
is derived from the Navier-Stokes equations (instead of from 
Euler's equations) given our choice of initial gas profile (see 
Appendix A). We neglect self-gravity and assume constant 
kinematic viscosity, v, and uniform temperature at all times 
(r « 10* K). We comment more on the astrophysic al m oti- 



vation behind our choice of initialization in section 5.1 



To obtain a constant u in the computation, the collision 
time r in the BGK scheme is set to ^^-^ (see discussion in 
section |t|) where U = 220 km/s, L = 20 kpc, Re^ « 750. 
The computation is performed on a 207 X 207 Cartesian grid 
and we follow the flow for about 12 Gyr. We would like to 
point out that in order to preserve the cylindrical symmetry 
of the flow on such a grid we had to use the high order limiter 
described in Ifuyn h (19 95). We specify boundary conditions 



in detail in section 5.1 



below. 



In figures ^ andg we plot the code's results for the time 
evolution of the gas surface density. Eg, the radial velocity, 
vr, and the ratio of the circular velocity at time t to its 



initial value. 



c{0) 



starting with an inviscid equilibrium so- 
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STEADY STATE ANALYTIC PREDICTIONS 
INITIAL CONDITIONS 

TIME EVOLUTION, CURVES SPACED EVERY 45,6 MYR 
PROFILE AT 2.28 Gyr 
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STEADY STATE ANALYTIC PREDICTIONS 
INITIAL CONDITIONS 

TIME EVOLUTION, CURVES SPACED EVERY 456 MYR 
EINAL PROEILE (ti>< 1 1 .84Gyr) 



R (kpc) 

Figure 3. Early viscous evolution (first 2.28 Gyr) of a disk of 
pure gas in a fixed gravitational potential. This time the initial 
radial velocity is set to the velocity profile given by the solution 
to the Navier-Stokes equation corresponding to the initial den- 
sity profile (see Appendix A for details). Coding for the different 
curves is given in the bottom panel of the figure. 



lution (i.e. with vr = Q initially )and a nearly uniform mass 
density distribution. In figures ^ and ^ we show a similar 
time evolution of these quantities but in this case we start 
from an initial vr consistent with the Navier-Stokes solu- 
tion corresponding to the initial density profile we choose *. 
Overplotted are the analytic predictions. These predictions 
are obtained by solving the Navier-Stokes equations for a 
compressible viscous fluid in cylindrical coordinates along 
with the equation of continuity: 



^3 [""R J. 



dVR 

dr 



dp 
dr 



dr 
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drj dvR 
dr dr 
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VR 
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Figure 4. Same run as in fig. ^ but with the results for the evolu- 
tion shown up to 11.84 Gyr. The steady-state analytic predictions 
(asterices) are compared to the BGK numerical results (thin solid 
curves and empty circles). 



MVR ^ vr \ 

\ dr r J 



drj / dvc vc^ 
dr \ dr r ) 
dT.a 



dr 



■Vr 



(6) 
(7) 



The fundamental assumption that permits the analytic cal- 
culation is that the circular velocity vc remains unchanged, 
and therefore one only needs to use the second Navier-Stokes 
equation (eq. ^) and the continuity equation (eq. ^). How- 
ever, this assumption can be tested a posteriori, when the 
solutions for density Eg and radial velocity vr have been 
computed. If one plugs these solutions in the first Navier- 
Stokes equation (eq. ™, it is easy to check that, to second 
order in viscosity it is satisfied. The analytic steady state 
solution of the system therefore is: 



vc{r) 
VR{r) 



r 



+ r^ 



2 I 2 

'■o+r 



(8) 
(9) 

(10) 



these solutions are given in appendix A 



One important remark is that this set of solutions de- 
pends only on the shape of vc through the parameter ro. 
Therefore, provided that the gravitational potential is fixed, 
the radial infiow of gas {vr is always negative) after a cer- 
tain time should be the same whatever the initial conditions 
we choose for the gas. Unfortunately, things are more com- 
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plicated for the density profile because tfie analytic solution 
for Sg(r) diverges for r = 0. However, there is an obvious 
comment that we can make: the steady state density profile 
is not exponential, which implies that it is the interplay be- 
tween viscosity and star formation that must be responsible 
for turning stellar profiles into exponentials in the viscous 
scenario. 

The first thing to notice in figures ^, |2| ^ and ^ is that 
the BGK scheme drives the gas towards a radial velocity pro- 
file which is in excellent agreement with the predicted ana- 
lytic profile, independent of the initial profile for vr. We will 
show in section 5.4 that this also holds when we start from 



a different initial density profile and/or a different gravita- 
tional potential profile, again in perfect agreement with the 
analytic expectations. The situation for the density profile 
is different because to obtain a good match to the analytic 



fit for a density profile which diverges like 



0, 



one requires an infinite amount of time and mass. As we 
can only run the computation for a finite amount of time 
and handle finite numbers, we obviously cannot achieve the 
steady state density distribution, even though one clearly 
sees that the gas will perpetually pile up in the center of 
the simulated disk. We also point out that the circular ve- 
locity of the gas, vc, remains unchanged up to a tenth of 
a percent throughout almost all of the disk, which numeri- 
cally confirms the validity of our analytic assumption that 
it is constant in time. Finally, although we do not show the 
results in this paper, we have also verified that decreasing 
Rcc, which is equivalent to increasing i/, results in the gas 
showing a similar radial fiow in shape, but with a higher 
normalisation. 

The success of BGK in giving viscous radial flows on 
the order of 1 km/s in a disk rotating differentially at 220 
km/s is remarkable. It is a technical success which is neces- 
sary for studying star formation on a viscous time scale in 
disks since chemical evolution models (Lacey and Fall 1985, 
Clarke 1989) find that radial flows on that order can explain 
the radial distribution of metals in disk galaxies. 



5 STAR FORMATION IN VISCOUS GALAXY 
DISKS 

5.1 Initialization and Boundary Conditions 

For our simulations with star formation, the initialization is 
similar to the initialization of the purely gaseous disk de- 
scribed in section ^ Here we elaborate on our initial and 
boundary conditions, and discuss our motivation for choos- 
ing this problem set-up. We initialize our gas to be in a flxed 
dark matter halo and we model the dark matter halo with 
the logarithmic potential defined by equation ^. We choose 
this form for $ because of its simplicity which allows us to 
perform analytical computations for comparison. Further- 
more this potential has the property that the rotation curve 
corresponding to it becomes asymptotically flat with veloc- 
ity Ho at large radii as is observed for many spiral galaxies. 
All the results presented in this paper are from runs with vq 
— 220 km/s. The total mass of the dark halo is chosen to be 
10^^ M0 and that of the gaseous disk 5% of the dark matter, 
i.e. 5 X lO^'^M©. Given that the gas is a small fraction of the 
dark matter, we neglect its self-gravity and assume that it 
evolves in the fixed potential created by the dark matter. 



We perform our simulations with an isothermal equa- 
tion of state, meaning that we assume that at the end of each 
iteration the gas instantanesouly cools to ~ 10*K. We choose 
this temperature because it is the lowest achievable from HI 
line cooling alone. Below that temperature we would have to 
worry about the formation of molecular hydrogen, which we 
neglect because we do not have the necessary resolution to 
follow the formation and fragmentation of molecular clouds 
in the context of an entire galaxy. Accounting for a molec- 
ular component would also require us to follow an atomic 
and molecular fluid simultaneously, which would consider- 
ably increase the complexity and computational cost of the 
simulation and therefore is beyond the scope of this paper. 
For a 7 = I gas a temperature of ~ lO^K corresponds to a 
sound speed of ~ 10 km/s. 

Since we want to initialize the gas with self-consistent 
viscous initial conditions, we set the radial flow to be given 
by the proflle derived from the Navier-Stokes equations (see 
Appendix A), assuming centrifugal equilibrium for a partic- 
ular initial gas density profile in a potential $ with a core 
radius ro (eq. Indeed, due to the observed thinness of 
stellar disks in nature (scale height ~ 200 pc), and the dif- 
ficulty of dissipating stellar motions perpendicular to the 
galactic plane, it seems reasonable to assume as an initial 
condition for the radial velocity that the gas collapsing to a 
disk in the dark matter halo was in approximate centrifugal 
equilibrium before it started forming stars. 

We choose 20 kpc for the disk radius of our simulations 
because based on observations of OB stars (Chromey 1978) 
and of HII regions (Fich and Blitz 1984), the stellar disk 
of our own Galaxy seems to have an outer cutoff around 
this value. Our computations are performed on a 207 by 207 
evenly spaced Cartesian grid, so that for a disk with radius 
20 kpc the cell size is about 200 pc per side. Again, we 
emphasize that we are not trying to resolve star formation 
on a small scale. We are trying to model the first order 
effect of turbulence, a phenomenon which might govern star 
formation on a large scale. We run each simulation for about 
12 Gyr. On a grid of 207 X 207 this takes approximately 
50,000 iterations and requires about 11 hours of CPU time 
per processor on an Origin 2000 with four rl2k processors 
working in parallel. 

We investigate the case where the kinematic viscosity, 
1/ is constant because it is a simple case for which we can 
compute analytic solutions for the radial velocity exactly 
and compare them to the results from our hydrodynamical 
simulation. As already described in section W, to achieve a 
constant u, we set the collision time r = . We again 

take U = 220 km/s, L = 20 kpc, Rec = 750 "and X = jpp 
with T ~ lO'^K. We underscore that our choice for u is 
motivated by our intention to model a disk with a critical 
Reynolds number ~ 750, because in this way as discussed in 
section ^, we model the lowest order effect turbulence has 
on a flow, namely to turn it into a more viscous fluid with 
Re — Rec . 

Our boundary conditions are as follows. Outside of a ra- 
dius of -Rdisk on our 2D Cartesian grid, we keep two "rings" 
of ghost cells, each "ring" being one cell thick. Beyond these 
ghost cells we do not follow the evolution of the gas. In this 
way we have carved a circular grid out of the square carte- 
sian grid. At the end of each updating time step of the code, 
we replace the values of the velocities in the ghost cells with 
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Figure 5. Comparison of the final stellar profiles (after f» 12 Gyr 
of evolution) for different values of the star formation timcscale 
to the viscous timescale, (3. The initial radial velocity profile is 
derived from the Navier-Stokes equation for the corresponding 
initial density profile, = Eo/(l + (^)^) where Vp = 78 kpc, 
So = 41.08Mq/pc^, and ro = 2.56 kpc in the gravitational po- 
tential (see Appendix A for details). 



the values of the closest active ring of cells (these for which 
the evolution is computed). This means that we are actu- 
ally performing a constant radial extrapolation to reset the 
values of all the macroscopic quantities (density, momen- 
tum and energy) in the ghost cells. Basically, this is an im- 
plementation of "outflow" boundary conditions, where the 
active region dictates the behaviour of the flow inside the 
boundaries. Nevertheless, despite this a priori free bound- 
ary condition, we find that the mass fiow across the outer 
boundaries of the grid is under control. At most our bound- 
ary conditions lead to a maximal loss/gain in total mass 
on the grid of 1.7% (see tables ^ and S). This was to be 
expected, since analytic calculations (eq^ and fig. ^ show 
that for the critical Reynolds number we are using, the outer 
radial fiow should eventually settle around a value smaller 
than 0.5 km/s. In the future, it will be interesting to inves- 
tigate the effect of infall onto the disk and for this we will 
obviously have to alter these boundary conditions. 



5.2 The star formation timescale 

The first point that we wish to verify is that if star formation 
proceeds on the same timescale as the viscous redistribution 
of mass and angular momentum in disk galaxies then the 
stars attain an exponential density profile. We therefore de- 
fine the viscous timescale, tvisc, by r/vR. This estimates the 
time it takes a disk annulus to move a radial distance r. Our 
star formation law then simply reads: 



(11) 



where, unless explicitly specified, t-i, = fvisc. In the steady 
state, tvisc ~ {ro + r^)/'^, so that at large radii the as- 
sumption of a constant kinematic viscosity corresponds to 
a star formation rate per unit area of Esfr oc Sg/r^. This 
SsFR corresponds to a Schmidt law with n = 1 modulated 
by r~'^ which means that the star formation rate at large 
radii is very slow. Note that since vr varies with space and 
time until the steady state is achieved, after which point it 
only varies with space, fvisc (and hence t*) has a similar be- 
haviour. We assume, as do Lin & Pringle (1987), that once 
the stars form, they are frozen out of the viscous evolution 
and simply move on circular orbits with a velocity given 
by the rotation curve at the location at which they form. 
Hence we neglect any process that might give the stars a 
velocity dispersion after they have formed. As a result of 
these assumptions, the implementation of star formation in 
the simulation only introduces a sink term, —T,g{r,t)/ti,, in 
the continuity equation ^ 

We perform a set of three runs with varying star for- 
mation timescales. Using Lin and Pringle's notation of /3 = 
t-k/tvisc, we run the code for the same three cases that these 
authors considered: (3 — .1, 1, 10. The initial gas density is 
= Eo/(l + (;^-)^) where Vp = 78 kpc, Eo = 41.O8M0/pc2, 
and ro — 2.56 kpc in the gravitational potential. Results are 
shown in figures ^ 6 and ^ 

Figures |^ and 3 show a comparison of the stellar pro- 
files we obtain after about 12 Gyr for these three different 
values of /3. Just as Lin and Pringle claimed, we find that 
the initial gas profile is nearly frozen into the stellar profile 
for small /3 because star formation is so rapid that almost 
no viscous redistribution of matter can take place before the 
stars form. Furthermore, as illustrated in the left panels of 
figure 1^ which show the time evolution of the gas density, 
radial velocity, and stellar density (top to bottom), almost 
all the gas is converted into stars by 12 Gyr (more than 
98 %). This is inconsistent with present-day spirals like the 
Milky- Way which still have a supply of gas of order 15 % 
(Prantzos & Aubert, 1995) to fuel ongoing star formation in 
their disks. 

At the other extreme is the case of large /3. Here too 
much viscous evolution occurs and the resulting stellar pro- 
file is more a power law than an exponential (see bottom 
panel of figure ^. Furthermore, the amount of stars formed 
during 12 Gyr is negligible (a few %), again in confiict with 
observations for Milky- Way type spirals which should have 
formed between 80 and 90 % of their stars. 

Only the l3 — 1 case, i.e. the case where stars form on 
the viscous timescale, gives both approximately an exponen- 
tial stellar profile (see residuals in top panel of figure ^ and 
a stellar density which is in better agreement with Milky 
Way values. The stellar density is still low for this simula- 
tion with only about 40 % of stars formed, but this number 
depends sensitively on the steepness of the initial gas profile 
(see table ^ and |^ . If we use a reasonably steeper (and prob- 
ably more realistic) initial profile we can bring the amount of 
stars formed in perfect agreement with the observations, e.g. 
for the case where E* = Eo/(l -I- (7-)''), rp = 5.12 kpc, and 
ro = 1.92 or 2.56 kpc, about 80 % of th e di sk's gas is turned 
into stars after « 12 Gyr (see section 5.4). Looking at the 
lowest middle panel of figure 0, we also see that after about 
6 Gyr (9th thin solid curve in the bottom middle panel, the 
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Figure 7. Time evolution of gas density (uppermost panels), radial velocity (second panels from top), tangential velocity (third panels 
from top) and stellar density (bottom panels) for /3 = .1 (left panels), i.e. a simulation where star formation takes place quickly compared 
to the viscous timescale, /3 = 1 (middle panels), i.e. a simulation where star formation occurs on a viscous timescale and /3 = 10 (right 
panels), i.e. a simulation where star formation is much slower than the viscous timescale. 



first one being the lowest), the final exponential profile is 
more or less established in the sense that there is little evo- 
lution in its scale length thereafter: the regularly spaced (in 
time) curves are much closer to one another than they are 
at earlier times, and their slope is not changing noticeably. 
The bottom panels of figure M lead to the general statement 
that the evolution of the stellar profile proceeds much more 
quickly at early times than at later times because the gas 
content is higher and the radial fiows are generally faster in 
the beginning of the simulation. 

In their models of star formation in a viscous disk Fir- 
mani, Hernandez & Gallagher (1996) also claim that there 
appear to be two distinct phases of star formation, the first 
corresponding to a stage of significant viscous mass and an- 
gular momentum redistribution, and the second correspond- 
ing to a stage where star formation has reduced the amount 



of gas in the disk to the point that viscous evolution is very 
slow and star formation proceeds slowly with an unchanging 
radial velocity profile for the gas. We remark that interpret- 
ing this as two distinct phases of star formation must be 
done cautiously because Firmani et al. (1996) initially take 
their gaseous disk to be in inviscid centrifugal equilibrium, so 
that when they begin to follow the disk evolution with equa- 
tions of motion which include viscosity, the gas first tries to 
find a (non-existent) equilibrium for the viscous case and 
this causes a lot of initial mass and angular momentum re- 
distribution. For instance, comparing our figures ^ and ^ 
we clearly see that the steady state for the radial velocity 
is reached much earlier when we start from an initial radial 
velocity in agreement with the Navier-Stokes solution corre- 
sponding to the initial density profile (middle panel of fig ^ 
than when we start from the inviscid equilibrium solution 
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Figure 6. Fitting residuals for /3 = 1 case (top) and /3 = 10 case 
(bottom). Thin solid lines are the residuals for the best power law 
fit to the stellar profile and the long dashed lines are the residuals 
for the best exponential fit. The fit is done only over the radii 
enclosed within the vertical short-dashed lines. The scale length 
of the exponential fit, r^, and the exponent of the power-law fit, 
Op, are given on each panel. 



(middle panel of fig ^). We also note that, whatever the ini- 
tial conditions, the time taken to reach the steady state is 
quite short (~ 2 Gyr for the longest initial inviscid equilib- 
rium case) compared to the duration of the whole simulation 
(about 12 Gyr). 

Since there is yet no absolutely convincing explana- 
tion for an equality between the viscous and star forma- 
tion timescales, it is worthwhile to explore how much leeway 
there is in the equivalence. Saio & Yoshii (1990) explored the 
cases where /3 is .5 and 2, as opposed to an order of mag- 
nitude greater or smaller than fi =\. They find the /3 = .5 
case to yield a profile close to the pure exponential disk 
(van der Kruit 1987) and the /3 = 2 case to yield a disk 
with a pronounced bulge and disk. They concluded that the 



star formation timescale does not have to exactly equal the 
viscous timescale to produce an exponential stellar profile. 
Although we do not show results in this paper, we have run 
simulations with different fi and confirm that provided it 
does not differ from unity by more than a factor 2, we get 
exponential profiles of a similar quality as those presented 
here. 



5.3 Dependence on initial gas profile 

The next step is to check the robustness of the t* ~ t^isc hy- 
pothesis to varying initial conditions. With this intention, we 
perform numerical experiments, keeping the rotation curve 
fixed (uo = 220 km/s and ro = 2.56 kpc) while varying the 
initial gas profile. We have tried the following two profiles: 
E^(r) = Eoe"^ and E^(r) = Eo/(l (;^)"'') where Eq is a 
constant, and are the scale lengths of the exponential 
and power law profile respectively and is either 2 or 4. For 
each of the runs with rotation curve prescribed by = 220 
km/s and ro — 2.56 kpc. Table |l| presents the values for Eq, 
rp, re the stellar exponential scale length, Va, resulting from 
an exponential fit to the stars after about 12 Gyr, the initial 
half mass radius, th, (i.e. the radius containing half of the 
mass in the disk) the exponent from a power law (PL) fit to 
the stellar profile, the percentage of initial gas that has been 
turned into stars after 12 Gyr, the percentage of initial gas 
that is left after 12 Gyr, the percentage of initial gas mass 
that has been lost from the grid or gained after 12 Gyr, the 
stellar density and the final gas density at the solar galacto- 
centric radius (8.5 kpc, e.g. Binney & Tremaine, 1987) after 
12 Gyr, and finally the initial gas density at 8.5kpc for com- 
parison to the sum of the stellar and gaseous densities at 
that location after 12 Gyr. 

We discussed the ru n co rresponding to the first entry 
of the table |l| in section 5.2. Figure ^ shows a comparison 
of the stellar and gas profiles for the other 4 runs presented 
in this table. We have overplotted the exponential fit to the 
final stellar profiles. We confirm that if the distribution of 
gas in a disk is originally exponential, and star formation 
proceeds on a viscous timescale, the exponential profile will 
be transfered to the stars. However due to viscous flows, the 
exponential stellar profile will have a different scale length 
than the scale length of the original gaseous disk. The final 
exponential scale length is about .5 th which is consistent 
with the expectation of Lin and Pringle's (1987) model. 

For stellar profiles which did not originate from a 
gaseous exponential profile, it is clear that the divergence 
from the exponential is large in the central regions where all 
the gas has been consumed rapidly, resulting in the initial 
gas profile being almost frozen into the final stellar profile. 
The central regions of real galaxies are also places where 
we certainly cannot neglect self-gravity and the presence of 
bulges and/or bars, which will have a dramatic impact on 
the gas. Therefore we exclude the inner 3-4 kpc from the fits 
of our exponential profiles. Additionally, to be conservative 
and minimize a possible contamination from the boundaries, 
we also exclude at least the outermost 2 kpc from our fits. 
The residuals, plotted in the right panel of Figure ^ clearly 
argue in favor of an exponential stellar profile over a power 
law one in most cases, although the departures from an ex- 
ponential are not negligible and can easily reach 20 % and 
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Figure 8. Comparison of stellar and gas profiles after 12 Gyr for varying initial gas profiles. From top to bottom the initial gas profile 
is = Soe"^, = So/(l + i:^)'^), ^ = So/(l + i^)'^), and EJ, = So/(l + C^)*). The values of the parameters Sq, r^, and 
are given in table |l| The parameters defining the gravitational potential are the same for each of these runs, namely ro = 2.56 kpc and 
vo = 220 km/s. 
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-2.41 
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25.1 
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PL, ap = 
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10.2 


3.47 
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-2.85 


59.6 


41.5 


-Hl.l 


45.2 


31.8 
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Table 1. Results from runs all having the same rotation curve {vq = 220 km/s and = 2.56 kpc) but different initial gas profiles (PL 
stands for power law with exponent ctp, and EXP stands for exponential. See text for more details). Column PL fitcontains exponents 
from the best power law fit to the final stellar profile; the following column gives the percentage of initial gas that has been turned into 
stars after 12 Gyr. From left to right, the next columns give values for the percentage of gas left after 12 Gyr, the percentage of mass 
lost (or gained) from the grid after 12 Gyr, the stellar density at the solar radius after 12 Gyr, the final gas density at the solar radius, 
and the initial gas density at the solar radius. 



more in some parts of the disk. Nevertheless, the profiles 
are generally exponential to within 10% over several scale 
lengths. Since our simulation grid only extends to 20 kpc, 
of which we only take a range spanning about 14 kpc seri- 
ously, it means that for stellar profiles with a scale length 
of around 3.5 kpc we only have a grid big enough to tell if 
the exponential fit holds over 4 scale lengths, which is not 
a significant improvement over Gunn's (1982) exponential- 
like profile. However for cases where ~ 2 kpc, we obtain 
exponentials which hold over more than 5 scale lengths (see 
table I and fig. 

On another note, in agreement with Lin & Pringle 
(1987) we find that the final stellar exponential scale length, 
rs , can be roughly approximated by half the initial half-mass 
radius, rn, for the gas (fig. ^). Hence, even though the fi- 
nal stellar profile is roughly exponential and in its shape 
forgets the initial gas profile, the magnitude of the stellar 
exponential scale length preserves the memory of the degree 
of central concentration of the initial gas distribution. 

After 12 Gyr we find that there is about 6 — 60% gas 
left in the galaxy depending on the steepness of the initial 
gas profile (see tables |l] and and obtain agreement with 
the final gas fraction of present day disks (~ 15 %) only 
for quite steep profiles. This is of considerable importance, 
because it means that by simply equating the star formation 
timescale to the viscous timescale, for disks with initially 
steep profiles we not only recover the exponential stellar 
profile, but we naturally get the star formation efficiency 
correct. We also find that our stellar densities at the solar 
galactocentric radius, Rq = 8.5 kpc are on the low side of the 
Milky Way value (35±5M0/pc^, Kuijken and Gilmore 1989) 
for steep profiles, but not completely incompatible with it, 
provided a moderate amount of infall takes place during the 
12 Gyr of evolution to boost the star formation rate. To our 
knowledge, ours is the first investigation that does not evolve 
the viscous disk in arbitrary time units, and then associate 
the present time with the time corresponding to the case 
where the stellar profile has solar neighborhood values at the 
solar radius (see for example Olivier et al. 1991, Ferguson 
and Clarke 2001). We stress that it is mandatory to have 
real physical units to be able to comment, as we do above, 
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Figure 9. The scale length of the stellar profile after « 12 Gyr is 
plotted against the initial half mass radius, rjj , of the gas profile. 
Asterices correspond to runs with core radius ro = 2.56 kpc for 
the gravitational potential <I>, but different initial gas profiles, and 
open trian gles correspond to runs with the same initial gas profile 
(cf. section 5.4) but a range of ro bracketted by 1.92 and 5.12 kpc. 



on whether the viscous timescale is a natural timescale for 
star formation. 
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Figure 10. Comparison of tlie stellar (left), gas and radial velocity (middle) profiles after « 12 Gyr for different rotation curves obtained 
by varying the core radius ro of the gravitational potential $ (from top to bottom r'o = 1.92, 2.56, 5.12 kpc). The initial gas profile for 
each of these runs is = Eo/(l + (7-)'*) with rp = 5.12 kpc and So = 402Mq/pc^. The residuals for the best exponential (dashed 
line) and the best power law (thin solid line) fits, within the radii delineated by the vertical short dashed lines, are given in the right 
hand panel. 



5.4 Dependence on rotation curve 

Lin and Pringle (1987) also stipulate that a stellar expo- 
nential profile always arises if t-^ « t^is^ regardless of the 
shape of the dark matter potential. By varying ro in the 
expression for the dark halo potential, $, we explore star 
formation for viscously evolving disks in potentials ranging 
from those whose rotation curves rise steeply in the central 
regions to those with more slowly rising curves. For each of 
the runs for which we present results, the initial gas den- 
sity profile is = Eo/(l + {^f) with = 5.12 kpc and 

Eq = 4O2M0/pc^. We again start from the initial radial 
velocity profile consistently derived from the Navier- Stokes 
solution (see Appendix A). 

Figure |l^ shows a comparison of the stellar profiles after 
about 12 Gyr for the cases where ro = 1.92, 2.56, 5.12 kpc. 
This figure also shows the corresponding radial velocity and 
gas profiles at about 12 Gyr. As mentioned in section M, we 
emphasize that the steady state radial velocity profile is now 



different for these three simulations, because it depends on 
ro. Once again, we draw the attention of the reader to the 
fact that the BGK scheme converges to the analytic solution 
for the radial velocity (even if less rapidly in the case of 
the most slowly rising rotation curve). Like Lin & Pringle 
(1987) who state that the results of viscous evolution are 
insensitive to the central concentration of the dark matter, 
we find that the stellar densities reach only slightly higher 
values and have slightly steeper scale lengths (rs = 2.07 kpc 
for ro = 1.92 kpc and 2.12 kpc for ro = 5.12 kpc) when 
the rotation curve rises more steeply. Once again, we are 
aware that for the high central densities we have here, self- 
gravity becomes important and we plan to tackle this issue 
in a subsequent paper. In this work we simply exclude the 
central region of the disk from the fit. 

Table ^ summarizes the findings of these runs. Even 
though now not all the rotation curves are good represen- 
tations of the Milky Way rotation curve, we take the solar 
galactocentric radius value of density as a rough estimate of 
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Initial Parameters 



Profile So ro re Ts th PL stars gas mass E» at Rq Eg at Rq Sg at Rq 

(M0/pc2) (kpc) (kpc) (kpc) (kpc) fit loss (Mo/pc2) (M0/pc2) (Mq/pc^) 



PL, Op = 4 402. 1.92 - 2.08 5.07 -3.78 82.4 17.6 -.02 28.2 17.4 45.4 

PL, Qp = 4 402. 2.56 - 2.08 5.07 -3.78 81.3 18.7 -.01 27.6 17.9 45.4 

PL, ap = 4 402. 5.12 - 2.12 5.07 -3.72 71.8 28.2 -f.Ol 24.4 21.0 45.4 



Table 2. Results from runs starting with the same initial gas profile (Eg = Eo/(l + (7-)*): where rp = 5.12 kpc) but different rotation 
curves obtained by changing the value of ro. vq, the asymptotic velocity is 220 km/s for all these runs. Columns are the same as those 
of table ry except rp (third column) has been replaced by ro . 



what we should find. We again see that we match the stel- 
lar surface densities in the solar neighborhood reasonably 
well and the quantity of gas left over in the disk after « 12 
Gyr is consistent with the observational estimate of « 15 %. 
Furthemore, the relationship between the stellar exponen- 
tial scale length and the initial half mass radius of the gas, 
{va = .5 r£f) still holds (see triangles in fig. |^. We also learn 
from figure y that for a fixed initial gas profile, the steepness 
of the rotation curve (quantified by ro in $) introduces very 
Httle spread in the correlation between initial gas half mass 
radius and the stellar exponential scale length after 12 Gyr. 

Another thing we can gauge from these runs, is how 
large an effect the magnitude of the radial velocity has on 
the final stellar density. Fig. |l^ and table |^ clearly show that 
the magnitude of the radial velocity is correlated to the total 
amount of stars formed at 12 Gyr. For example the run with 
the most slowly rising rotation curve (ro = 5.12 kpc) has a 
lower steady-state radial velocity throughout the disk than 
the runs with steeper rotation curves. In the steady state, 
the maximum value of vn for the run with ro = 5.12 kpc 
is a factor 1.8 smaller than the value for the ro = 2.56 kpc 
run. Since ti, ~ tvisc, where tvisc = t/vr, the star formation 
timescale is longer for the run with ro = 5.12 kpc and hence 
at 12 Gyr it has formed 10 % fewer stars than the run with 
ro = 2.56 kpc. This points to another important consider- 
ation, namely the value we take for the critical Reynold's 
number, Rec, which determines the amount of viscosity and 
hence the magnitude of the radial fiow. We took Rcc ~ 750, 
but as we discussed in section ^ it could range from 10— 10'^. 
If we had taken Rec ~ 100, for example, then the collision 
time in the scheme, and hence the kinematic viscosity would 
increase by a factor 7.5 causing the radial infiow velocities to 
have a maximum velocity of about 7.5 km/s. Given the lack 
of observations which pinpoint the radial velocity of gas in 
slowly evolving disks, we rely on arguments concerning ra- 
dial metallicity distributions (Lacey and Fall (1985), Clarke 
(1989)) to constrain the magnitude of the radial velocity to 
about 1 km/s, but we underline that this is indirect evidence 
for the magnitude of vr and hence it is still debatable. 

Finally we remark on the gas profiles. Consistent with 
Lacey and Fall's (1985) results, we find that the gas surface 
density profiles at about 12 Gyr are deficient in gas at radii 
less than about 5 kpc (Figs. g| and |l^). Also in general we 
find that the gas profiles after 12 Gyr for radii greater than 5 
kpc, if fit by an exponential, have a larger scale length than 
the stellar profile. This is also consistent with Lacey and 



Fall's results. Nevertheless, we do not want to overinterpret 
our results for the gas since the gas profiles can be strongly 
infiuenced both by feedback effects induced by supernovae 
or stellar winds and infall. 

5.5 Star Formation Rates 

Another test of the viscous model for star formation in galac- 
tic disks, is the comparison of the star formation rate with 
the gas surface density as a function of time. 

On figure we plot the results from our model with 
an initial gas profile given by Eg — So/(l + (7-)^), with 
rp = 5.12 kpc, and evolving in a potential characterized by 
ro = 2.56 kpc and no = 220 km/s (other results for this 
profile were shown in fig. |l^. On the same plot we also in- 
dicate the observational constraints set by Kennicutt's Ha 
data (Kennicutt 1998) for the SFR surface density Espr and 
CO data for the gas surface density, Egas- Top panel of fig. ^ 
gives results of the early evolution (up to 1.2 Gyr), whereas 
bottom panel of the same figure shows the late evolution 
(from 1.6 to 12 Gyr). There are two important remarks to 
make concerning this figure: the first is that one can clearly 
see star formation happening inside out, with the central 
region having the highest star formation rate in the begin- 
ning. The latter progressively decreases due to the lack of gas 
to fuel it at late times. But this monotonic behaviour only 
takes place after the radial velocity profile has settled down 
to the steady state solution, because we see in the top panel 
of fig. ^ that there are some regions with a high density of 
gas that form a low amount of stars in the early stage. This 
is due to outfiows of gas for which we quench star formation. 
We explicitly only allow stars to form where the radial ve- 
locity is directed towards the center of the disk. The second 
point is that although we naturally reproduce the slope of 
the relationship and the cut-off at low gas surface densities 
seen in the data, the locus occupied by our curves at late 
stages is offset towards low values. This, however, should not 
be considered as a failure of the model, as the normalisation 
of the total gas mass is somewhat arbitrary. Had we started 
the simulation with 15 % of baryons in our dark matter po- 
tential instead of 5 % (which is also a reasonable cosmolog- 
ical value for the baryon fraction), the star formation rates 
and gas surface densities would have increased by a factor 
of three, bringing the model in perfect agreement with the 
observational data. Also, the picture would be more com- 
plicated if we introduced even a moderate amount of infall. 
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mation has to proceed on a viscous timescale in order for 
galaxies to end up with exponential light profiles, so the re- 
distribution of angular momentum cannot be very efficient. 
Although in the simulations presented in this paper, we do 
not take into account feedback via supernovae explosions 
or stellar winds that would reheat and redistribute the gas, 
we show that viscous evolution does not lead to a dramatic 
decrease in the specific angular momentum of the gas. 

On figure ^ we plot the time evolution of the cumu- 
lative specific angular momentum, j'c = J /M, of both gas 
and stars. This quantity is obtained by computing the total 
amount of angular momentum in concentric gaseous/stellar 
disks and dividing it by the gaseous/stellar mass in these 
disks. There are a couple of things that we would like to 
emphasize on this figure. First is that in all the runs, the 
stellar angular momentum distribution closely follows the 
initial gas angular momentum distribution, which already 
points out that viscous redistribution of angular momentum 
must be small. This is particularly noticeable in the inner 
regions of the disk. Second, one can immediately see from 
figure ^ that the final gas angular momentum distribution 
is always higher than the initial one, especially at large radii. 
In fact, there is a clear trend that the steeper the initial gas 
density profile and rotation curve, the larger the discrep- 
ancy between the initial and final specific angular momenta 
of the gas. The reason for this is the following: since ini- 
tial gas density profiles are generally decreasing functions 
of radius, equation |Tl| combined with equation ^ imply that 
the star formation rate is much higher in the center of disk 
galaxies. Hence star formation will proceed inside out, con- 
suming gas with low specific angular momentum first. As 
a result, after 12 Gyr, the leftover gas is naturally located 
in the outer regions of the disk where the specific angular 
momentum is higher in general, and the viscous timescale 
so large that redistribution is minimal, naturally yielding an 
increase in the cumulative specific angular momentum. 
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6 CONCLUSIONS 



Figure 11. Time evolution of the star formation rate surface 
density as a function of the gas surface density in one of our 
simulations defined by = So/(l + (f-)*); with Vp = 5.12 
kpc, ro = 2.56 kpc, and vo = 220 km/s. Solid lines are our star 
formation rate profiles at times indicated on the figure. Shaded 
region is observational data by Kennicutt (1998). 



because both star formation rate and gas surface density 
would also go up. Furthermore, as discussed at the end of 
section oA, had we taken a higher value for the viscosity, 



we would have found faster star formation timescales which 
would also have increased the star formation rates. 



5.6 Angular momentum 

Finally, for the rotating viscous disk, one might worry that, 
as pointed out by Lynden-Bell and Pringle (1974), the equi- 
librium is reached when a portion of fluid with negligible 
mass carries all the angular momentum to infinity while the 
rest of the gas collapses to the center. This would only ex- 
acerbate the angular momentum problem seen in N-Body 
-f SPH simulations, if it were not for the fact that star for- 



We have performed the first numerical hydrodynamical sim- 
ulations probing the viscous evolution of isothermal and non 
self-gravitating disk galaxies embedded in static dark halo 
potentials. 

Despite the simplicity of such a model, where we have 
also neglected the return fraction of gas through stellar evo- 
lution, infall, and feedback, we obtain disks with approxi- 
mately exponential stellar profiles, regardless of their initial 
gas profile and rotation curve as long as ~ t-uisc- More- 
over, we obtain fractions of gas and stars after ~ 12 Gyr 
of viscous evolution that are consistent with observations, 
without having to invoke an efficiency parameter for star 
formation. Although a detailed explanation for the connec- 
tion between the viscous timescale and star formation is yet 
to be worked out in details, this strongly suggests that the 
viscous timescale might indeed be the natural timescale for 
star formation. It is also a more appealing alternative than 
fine-tuning the initial conditions of proto-galactic clouds des- 
tined to form disks in order to obtain exponential stellar 
density profiles. 

The specific results of our investigation are: 

• When t* ~ tvisc, the resulting stellar profile is approx- 
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Figure 12. Comparison of the cumulative specific angular momentum profiles for the initial gas distribution (solid lino) and both the 
final stellar (asterices) and gas (open circles) distributions, i.e. after 12 Gyr of viscous evolution. The dashed line shows the initial gas 
profiles. The exact analytic shapes of these profiles are given in tables |l| and |. The core radius, ro, of the gravitational potential is 
written in the lower right corner of each panel and vq, the asymptotic circular velocity is 220 km/s in all simulations. 
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imately exponential. When ti, is 10 times smaller than the 
viscous timescale then the initial gas profile is frozen into 
the stellar distribution and when t* is 10 times larger than 
the viscous timescale then the stellar distribution is better 
described by a power law. 

• Exponentials fit reasonably well the stellar profiles after 
~ 12 Gyr, regardless of the initial gas profile or the rotation 



curve. Hurtnermore, tnere is no requirement on conservation 
of specific angular niomauLum Lhi'oaghouL HLiii' formiiLlou. — 

• For rotation curve profiles similar to Milky Way profiles, 
(ro = 2.56 kpc and vo = 220 km/s in $), after ~ 12 Gyr 
the stellar densities in the solar neighborhood and the final 
gas fraction in the disk are consistent with observations, 
provided the initial gas density profile is quite steep. This 
suggests that the viscous timescale gives the star formation 
efficiency naturally. 

• The disk scale length, Vs, depends most strongly on 
the degree of central concentration of the initial gas density 
profile (figure The steepness of the dark halo rotation 
curve introduces very little spread in the initial gas half mass 
radius, rn, versus rs relation. 

• The cumulative specific angular momentum of the final 
distribution of gas is always higher than the initial cumu- 
lative specific angular momentum of the gas due to inside 
out formation. The steeper the initial gas profile, the more 
important is the discrepancy between the cumulative spe- 
cific angular momentum of the gas after ~ 12 Gyr and its 
initial cumulative specific angular momentum. The cumula- 
tive specific angular momentum of the stars after ~ 12 Gyr, 
on the other hand, closely matches the initial cumulative 
specific angular momentum of the gas. 

• The slope of the star formation rate density versus gas 
surface density in the viscous disk models is in fair agreement 
with the observed one, and the cut-off seen at low gas surface 
densities in the observations is naturally reproduced. 

Clearly, this work is a necessary first step which proves 
the ability of the BGK scheme to follow the physical pro- 
cesses occuring in a galactic disk. In following works, we 
shall focus on implementing more realistic conditions for the 
starting point of the simulation and the mechanisms driving 
the disk evolution. We also plan to extend the code to 3D 
and include self-gravity of both dark halo and gas. 
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APPENDIX A: INITIAL CONDITIONS 

Here we give analytic expressions for the initial radial ve- 
locities, Vfi, which correspond to the initial density profiles. 
Eg we use in the paper, i. e. knowing vc from equation H, we 
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compute u^j by solving the second Navier-Stokes equation 
(eq. 1). This yields: 

• if Sg(r) = Eo = constant, 

. ifE^(r) = Eo/[l + (;^f] 

''^^ ^^ {rl+r^){2rl + r^){rl + r^) ^^^^ 
. ifE^(r) = Eo/[l + (^)*] 

w^_ r(4r>g + rV^-3r^) 

(rg+r2)(2rg+r2)(r4 + r4) ^"^^^ 

. if E^(r) = Eoexp(-r/r,) 

^,^ r(4r,rg +r-V,^ -rrg -r^) 

Note that these initial conditions for the radial velocity 
have quite different asymptotic properties. For instance, the 
initial radial velocity corresponding to the constant density 
profile (eq is always negative, meaning the gas is con- 
stantly flowing in as the steady state radial velocity is also 
always negative (eq|^. However, both radial velocities go 
to zero like at large radii for the initial power profiles 
(eq and A3), which means that there should be a small 
initial outfiow of gas in these simulations. Finally the initial 
radial velocity prof ile c orresponding to an initial exponential 
density profile (eq A4 ) can have arbitrarily high outfiows at 
large radii, depending on the value of r^: the smaller the 
value, the more important the outfiow since oc rj^. 



